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A mesoscopic model for shear plasticity of amorphous materials in two dimensions is introduced, and studied 
through numerical simulations in order to elucidate the macroscopic (large scale) mechanical behavior. Plastic 
deformation is assumed to occur through a series of local reorganizations. Using a discretization of the mechanical 
fields on a discrete lattice, local reorganizations are modeled as local slip events. Local yield stresses are randomly 
Cq \ distributed in space and invariant in time. Each plastic slip event induces a Eshelby-like long-ranged elastic stress 

■^1- ■ redistribution. Focusing on quasistatic loadings and zero-temperature limit, extremal dynamics allows for recov- 

£S) ' ering many of the complex features of amorphous plasticity observed experimentally and in numerical atomistic 

simulations in the quasi-static regime. In particular, a quantitative picture of localization, and of the anisotropic 
strain correlation both in the initial transient regime, and in the steady state are provided. 



Nous presentons un modele numerique a l'echelle mesoscopique de la plasticite des amorphes. La deformation 
plastique est supposee resulter de series de rearrangements locaux de la structure. La discretisation sur reseau 
des champs mecaniques permet d'assimiler ces rearrangements a de simples glissements locaux. Les seuils de 
K^ . contraintes associes a ces evenements sont aleatoires et independants dans l'espace comme dans le temps. Chaque 

evenement plastique induit un champ de contrainte elastique a longue portee (type Eshelby). S'interessant au 
comportement quasistatique et a la limite d'une temperature nulle, l'utilisation d'une dynamique extremale per- 
met de rendre compte de nombreux phenomenes observes experimentalement et numeriquement dans le regime 
quasi-statique. En particulier nous presentons des resultats quantitatifs sur le comportement de localisation et 
l'anisotropie des correlations de la deformation plastique. 
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1. Introduction 

The prominent mechanical property of glasses is their brittleness. Nevertheless, the question of their 
possible plasticity has been discussed very early[l,2]. Limited to small scale or confined geometries (due 
to their brittleness), the plastic behavior of glasses is quite different from its crystalline counterpart. In 
particular, a notable pressure-dependence is often encountered and, in addition to shear flow, it is common 
to observe permanent densification[3,4]. Despite this early interest for the plastic behavior of oxide glasses 
and the primary importance of this phenomenon for mechanical contact properties, the subject has long 
remained poorly explored. The recent development of metallic glasses has given a new momentum to this 
subject [5]. Indeed, the mechanical strength of metallic glasses appears to be limited by shear-banding 
induced failure phenomena. 

One of the most important issues that constitute a barrier in the development of those materials is 
the proper mastering of size effects. Indeed, if the yield limit appears to be extremely high as compared 
to crystalline materials with the same chemical composition, yet failure appears to be quite brutal for 
large scale specimens, akin to brittle fracture. Amazingly, a closer analysis show that these materials 
may withstand a large number of macroscopic slips over different shear planes prior to failure. The proper 
understanding of the initiation of inhomogeneous plastic strain is thus essential to enhance the mechanical 
properties of these materials. 

In terms of theoretical modeling, a major difficulty long lied in the absence for amorphous materials of 
any obvious microscopic alternative concept to dislocations for crystalline plasticity. The understanding 
that plastic deformation resulted from a series of local structural reorganizations [8,9], commonly referred 
to as "shear transformation", (ST), was a major theoretical breakthrough introduced by Argon. This 
new paradigm, together with the fast development of computing facilities has opened the way to direct 
simulations of shear plasticity of glasses at atomistic scale (such as Molecular Dynamics) and mesoscopic 
scale (intermediate between the size of a structural reorganization and the continuum scale). 

The phenomenology of amorphous plasticity that emerged frome these numerical studies to be discussed 
in greater detail below (see also Ref. [10] for a review) is rather complex and questions the traditional 
view of plastic flow as smooth and regular. In particular phenomena such as intermittence, avalanche or 
quasi-localisation are shown to occur over a wide range of time- and length-scales. The observation of 
such scaling behaviors naturally calls for a description which necessarily includes fluctuations and not 
simply averages, even at macroscopic scales. 

In the framework of the complex phenomenology of amorphous plasticity, our ambition is thus the 
construction of the minimal model that reproduces a number of macroscopic features. By minimal, we 
mean that the least number of free parameters and computational complexity is our guide. For this we 
will have to sort between important features which survive at a macroscopic scale, and others (which 
may indeed exists) which will be disregarded if they do not impact the large scale behavior. Following an 
early work by Baret et al. [11], our model will be based on the competition between the respective effects 
of the structural disorder of the amorphous matter and of the elastic coupling induced by the reaction of 
the elastic bulk to local reorganizations[12[. 

In this paper, a particular emphasis is put on the geometrical features of the plastic strain field. We 
first detail the ambition and the limits of our modeling approach in section 2. The model is constructed 
in Section 3, the stress-strain relation is discussed in Section 4, the localization behavior is discussed in 
Section 5. Section 6 gives a quantitative characterization of the anisotropic strain correlation. A final 
discussion, Section 7, concludes this article. 
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2. Amorphous plasticity: from microscopic scale to macroscopic scale 

In this section we briefly discuss the complex phenomenology of amorphous plasticity that has emerged 
from the large body of numerical simulations performed in recent years and we introduce our modeling 
strategy. 

2.1. Molecular Dynamics 

In view of the obvious numerical limitations of atomistic simulations, most of the recently simulated sys- 
tems consisted of 2D Lennard-Jones model glasses or close variants [9,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28] 
which already capture most of the difficulties of the problem, and yet allowed to consider large system 
sizes. A few three dimensional works [29,30] indeed confirmed that little new physics emerged from the 
third dimension. Some simulations even made use of more "realistic" potentials for metallic glasses [31], 
amorphous silicon, [32,33,34,35], vitreous polymers [36,37] or silica glass [38,39] reaching similar conclu- 
sions. 

Although a wide consensus exists on the scenario of ST initially proposed by Argon [8] as validated by 
the numerous numerical studies [13,32,15,16,33,34,20,24,25,35,40,41], the quantitative characterization of 
these local reorganizations as well as their links to the local structure remains quite elusive. 

Beyond the question of the microscopic mechanisms, a special attention has been devoted to the com- 
plex character of the plastic behavior of glasses. At a microscopic scale, the "non-affine" corrections to the 
elastic behavior have been thoroughly discussed [13,15,16,17,19,20,38,42,24,25] as well as their possible 
consequences on the initiation of plastic reorganizations [25] . Strain localization has received a particular 
interest [14,18,30,26,43] as well as avalanche statistics [15,27]; akin to other non-linear mechanical behav- 
iors in random media (e.g. fracture), plastic deformation is not continuous but proceeds by successive 
bursts. 

2.2. Meso-models 

In parallel to the above discussed atomistic simulations, several works have developed numerical models 
at a mesoscopic scale, i.e. intermediate between the scale of an individual reorganization and the macro- 
scopic scale where continuum mechanics gives a faithful description. Assuming that plastic deformation 
resulted from a series of reorganizations, Argon and Bulatov early proposed a discrete model of amor- 
phous plasticity and studied its behavior at different temperatures [44,45,46]. Homer and Schuh [47] very 
recently proposed an extension of this model to follow the dynamics of plastic deformation of amorphous 
metals under shear using a kinetic Monte-Carlo algorithm. 

Along similar ideas, Baret et al. [11] developed a model of amorphous plasticity at zero temperature in 
antiplane geometry with two main ingredients: a structural disorder and a long-range elastic interaction 
to account for the stress redistribution due to local plastic reorganizations. This model, which exhibits 
critical properties, can also be considered as belonging to the same vein as other statistical modeling 
for earthquakes [48,49] or elastic line depinning [50,51]. In a close spirit, Picard et al discussed complex 
spatio-temporal behavior of a yield stress fluid [52]. In this model, the authors considered a quadrupolar 
elastic interaction but no other source of disorder than the one caused by the relaxation dynamics. 
Following similar ideas, Bocquet et al developed an analytic model for the elastoplastic dynamics of a 
jammed material taking into account the elastic stress redistribution[53]. Lemaitre and Caroli[54] recently 
discussed the avalanche behavior within a mean field model. In particular, they proposed to link ideas 
of effective temperature discussed by Sollich et al. [55] to describe the rheology of soft glasses and the 
internal stress fluctuations by varying the distribution of their elastic interaction. Jagla [56] discussed 



the effect of a local relaxation processes on the localization behavior. Sekimoto[57] discussed the internal 
stress as encoding past mechanical treatments. In a more traditional dislocation plasticity spirit, Zaiser 
and Moretti[58] discussed the avalanche behavior in a mesoscopic model very close to the above discussed 
ones. More recently Dahmen et al[59] proposed a variant of the mean field Ben-Zion-Rice earthquake 
model[49] incorporating a systematic weakening or hardening effect. 

2.3. Macroscopic scale: softening vs hardening 

Some of these studies have been dedicated to the characterization of the plastic behavior within a 
continuum mechanics framework [36,37,31]. 

The large scale behavior is observed to approach a so-called elastic-perfectly plastic constitutive law, 
where the macroscopic stress becomes independent from the macroscopic strain. However, the latter law 
lies at the edge of mechanical stability. Would stress increases with strain ("hardening" behavior) that the 
macroscopic strain field would tend to be homogeneous and deterministic at large scales (under a uniform 
stress) amenable to homogenization[60,61]. In contrast, a "softening" behavior (where stress decreases with 
strain) induces at large scales localization modes where strain concentrates on shear bands. In the latter 
case, homogenization is inappropriate, and the macroscopic behavior results from microscopic details 
(initiation and growth of the shear band) without much hope for a picture which would uniformly hold 
for a large class of amorphous materials. 

The perfectly plastic behavior lies at the borderline between these two classes of behavior. Hence, one 
cannot state any general results about the homogeneity of the strain field or the fate of a localization 
mode. It is necessary to resort to the scaling features of a meso scale modeling as it will determine the 
answer to such a question. Interestingly enough, this perfect plasticity character will justify the fact that 
fluctuations will survive up to the macroscopic limit, without scale cut-off. Hence, in contrast to many 
other classes of constitutive law, an enrichment capturing fluctuations at all scales is a need. 

2.4. Modeling strategy 

The vanishing of the tangent macroscopic stiffness is thus the key feature that allows for the relevance 
of more microscopic details than in hardening behaviors. However, not all details matter, and sorting 
out the important ingredients needed to enrich an effective macroscopic description is an important issue 
which has not been addressed as such in the abundant above mentioned literature. 

In conclusion, a clear view of up-scaling is still missing. This question is limiting inasmuch as the 
identification of a proper continuum description framework is not settled (e.g. is the plastic strain field 
uniformly distributed or rather localized in the thermodynamic (large size) limit?) and hence the question 
of what microscopic features survive at a macroscopic scale (and should be extracted from extensive 
microscopic studies) is not elucidated yet. The present work aims at addressing the question of the 
asymptotic homogeneity or heterogeneity of the plastic strain over various space and time scales. In order 
to approach this problem, we will consider a large scale modeling and argue for the relevance at large 
scales of the elementary mechanisms in line with microscopic or mesoscopic modeling. 

We will consider a large scale modeling and argue for the relevance at large scales of the elementary 
mechanisms in line with microscopic or mesoscopic modeling. Our model will be seen to reproduce the early 
work by Baret et al [11] in the context of anti-plane loading. The latter case corresponded to a parallel 
displacement field along say the z-axis, whose amplitude depends only on the transverse coordinates 
U z (x,y). We will here rather focus here on a different symmetry, e.g. plane strain, where U z = 0, and 
the in-plane displacement depends only on x and y. Because the modeling is proposed based on large 
scale asymptotics, we will retain only the most salient effects which have a chance to survive at a large 



scale and hence, our objective is to propose the minimal level of complexity. Nevertheless, a detailed 
comparison of scaling features with microscopic modeling such as that of Maloney and Robbins[26] can 
be performed and will be shown to be extremely close. Alternatively, such a mesoscopic model should 
not be used to address questions pertaining to the microscopic scale. Therefore, without questioning their 
validity, a number of microscopic properties reported in the abundant literature will not be retained if 
we do not anticipate a macroscopic manifestation of these observations. These choices should not be 
interpreted as a negation of a more complex microscopic reality than the one we propose here. On the 
contrary, our goal is to identify the most elementary model, based on some accepted observations, which 
is amenable to reproduce the large scale features of amorphous media plasticity. Part of this motivation is 
to be able to investigate asymptotic scaling features. For instance, the question of the asymptotic nature 
of the plastic strain distribution (homogeneously distributed as for hardening behaviors, localized as for 
softening behaviors, or a more subtle intermediate regime to be explicited) is at the very heart of our 
approach. 



3. Definition of the model 

The studied model is here constructed with the objective of capturing the large scale behavior of amor- 
phous media. In order to limit the computational cost of studying the model, only the two dimensional 
case of plane strain will be considered. It is recalled that similar behaviors where observed in molecular 
dynamics simulations in two and three dimensions, and hence this restriction is not limiting. As argued in 
the introduction, the macroscopic loading may induce changes in volume which are coupled with the acti- 
vation of plastic flow. However, for extended shear, no macroscopic dilation or shrinkage can be expected 
and hence the dominant macroscopic plastic strain will be purely deviatoric. Note that a macroscopic 
volume change may occur in the transient early deformation regime, but is not accounted for in our model. 
Because only the transient regime would be affected, the scaling analysis of the steady plastic flow should 
not be affected. Hence for simplicity pure shear boundary conditions are considered, with fixed principal 
axes in order to reach large strains. More precisely the mean total strain field is controlled as 




e)=e\ (1) 



In the following a discretization will have to be introduced to study numerically the model. A square 
lattice with bi-periodic boundary conditions will be used to avoid a systematic spatial bias induced by 
rigid boundaries. 

The microscopic structure of amorphous media is intrinsically heterogeneous, and hence, even in the 
absence of plastic reorganization, the elastic behavior of the material is heterogeneous. However, upon 
coarse-graining, this heterogeneity will vanish giving rise to an homogeneous behavior for which elastic 
constants can be derived from classical homogenization techniques. Thus, the elastic regime will be de- 
scribed as homogeneous. Moreover, the amorphous medium will be assumed to be isotropic. Although 
this has been questioned on the basis of molecular dynamics simulation for amorphous silica [39] , such 
an anisotropy is weak and does not affect the main features (and in particular scaling properties) of the 
model at large scales. Since no volume change will be considered, only one elastic constant, the shear 
modulus, /i, has to be introduced, which in a dimensionless model can be set to 2/i — 1. Note that this 
statement does not contradict the existence of elastic heterogeneities, and their role in the triggering of 
ST, but it merely says that those heterogeneities do not survive at a large scale. Thus we may simply 
resort to a deterministic homogeneous description at a large scale. Hence whenever we refer to a local 



stress, it will refer to that of the homogenized picture. The actual local stress would be related to the 
homogeneous one with a localization operator which is not introduced explicitly because it is not needed 
in the sequel. 

The shear transformations are the only source of non-linearity in the model, and hence at any stage, 
the medium can be unloaded and reloaded by the mere superposition of a uniform shear stress field (with 
a fixed orientation). Hence, the tensor stress field can be characterized by its sole deviatoric norm (or 
equivalent von Mises stress). Scalar (equivalent) stress a = a xx — o~ yy and strain e = e xx — e yy will be 
used (although we use a plane strain assumption, the conditions specified below will guaranty that plane 
stress also hold). The latter scalar stress component will be denoted as a, and called "stress" for simplicity. 
Similarly a single scalar strain component will matter, and will be denoted as the strain, e, in the sequel. 
With those notations the incremental elastic law reduces to the elementary form Aer = Ae. 

As above discussed, in amorphous materials such as glasses, pastes or foams, a common assumption 
consists in describing the macroscopic deformation as deriving from a succession of localized reorganiza- 
tions at some microscopic scale[9]. These regions can be seen as the coarse grained Shear Transition Zones 
(STZ) [9] whose sizes reach at most a few tens of atoms [21] in molecular dynamics simulations. The de- 
tails of such local rearrangements depend obviously on the precise structure of the material under study. 
They involve local change in the topology of atoms, grains or foam cells, and significant non-linearities 
either of geometrical or constitutive origin take place. The important feature is that these reorganizations 
are local is space. We discuss below the criterion to initiate such an event (Sect. 3.1) and its effect on the 
stress field (Sect. 3.2). 



3.1. Onset of local plastic transformation 

In order to describe at which stage of loading such a shear transformation will take place, a criterion has 
to be proposed. It is natural in the present framework to characterize locally the onset of a transformation 
through the local stress. Here local refers to the homogenized description. In addition, it is also at the 
discretization scale. However, our discretization is itself at a mesoscopic scale not to be confused with the 
microscopic one. Let us only note here that the full characterization of the local reorganization occurring 
under shear in amorphous materials is still an ongoing issue. In particular the existence for amorphous 
materials of a local stress threshold at atomic scale has recently been debated [24,62]. While the definition 
of a satisfactory criterion can be discussed at atomic scale, there is no doubt however that under sufficient 
coarse graining, one should recover a criterion based upon the homogenized stress tensor whatever the 
localization tensor. Irrespectively of the precise characterization of local instabilities at atomic scale, we 
assume in the present mesoscale model, that the coarse-graining is performed at a large enough length 
scale to allow us to safely use local yield stress criteria. The criterion for yielding characterizes the local 
configuration of atoms, and hence will display some variability. A local yield threshold for each discrete 
site x as a 1 {x) is introduced, and will be treated as a random variable in the sequel. For all sites, the 
same statistical distribution will be used, chosen for simplicity as a uniform distribution over the interval 
[0, 1]. Other distributions could be considered but are not expected to alter the generic behavior of the 
model beyond a few mesh sizes (and this has been checked numerically). 

It is to be noted that since no stress scale has been introduced in the model so far, rescaling the 
distribution to [0,?i] will only multiply all stresses by the same factor ?i. Similarly, the interval may be 
translated to [?Oi?o + £1], by adding a constant <^o to all stresses. Thus the linearity of the model can be 
exploited to match any yield stress interval. More importantly, other distributions, such as a Gaussian, 
behave similarly, and at a large enough scale, the local yield distribution can no longer be read from 
microscopic observables, as a consequence of universality. 

6 



3.2. Effect of a local transformation 

Beyond the confined region where the rearrangement takes place, a stress perturbation has to be ac- 
commodated by elastic strain throughout the material. This perturbation induces internal elastic stresses 
a e i which vanish with the distance, as dictated by linear elasticity (of a homogeneous medium as above 
argued). The general solution of elastic field induced by a localized perturbation can be written easily 
as discrete infinite series (multipolar expansion [63]). An important feature of those modes is the ab- 
sence of length scale. Therefore, these stress influence functions can be extrapolated down to a point-like 
singularity. As such, they are genuine Green functions attached to the elementary rearrangements. The 
constraint of not imposing externally any force or torque selects fields which decreases with the distance 
to the transformation zone as r~ n where n > D where D is the space dimension. The leading perturba- 
tion term is hence given only with the n = D terms. Three such modes exist in two dimensions, D = 2, 
which coincide with the external problem of a plastic strain imposed uniformly in a circular inclusion, the 
so-called Eshelby problem [12]. One is a pure dilation which corresponds to the expansion or shrinkage of 
the transformation zone. The two others are pure shear modes with different principal axes. Because the 
macroscopic loading is a pure shear, and large plastic strains are considered, the dilation mode will vanish 
on average. Thus, any change of volume will be discarded. Similarly, on average, the orientation of the 
deviatoric mode will be aligned with the macroscopic loading, and fluctuations around this orientation 
will be neglected. Consequently, only one mode will be retained to describe the elastic effect of a localized 
transformation. Let us underline that neglecting the two other Eshelby modes is a simplifying assump- 
tion. Introducing those modes would require a full tensorial framework for stress and strains, and hence a 
much more complex description for triggering an elementary plastic event. Mesoscopic scale models such 
as [11,52] thus consist in studying the effect of these dominant terms of the elastic interaction on the dy- 
namics of localized plastic events, independently of the microscopic details. It is stressed that there is no 
upper limit on the discretization scale, since the considered external Eshelby mode does not involve any 
specific length scale. However, if the discretization scale changes, the corresponding statistical features of 
the local rearrangements have to be coarsened at the appropriate scale and hence the microscopic disorder 
has to be adjusted to match a specific material. This question will not be addressed in the present study. 

An elementary local transformation zone occurring at point Xq is thus described by a stress field 
rjG(x — Xq) whose scalar amplitude r\ is the only relevant feature. Note that the translational invariance 
which is natural in an infinite medium still holds for the chosen bi-periodic boundary conditions in a 
finite size system. The precise form of G is discussed below in Sect. 3.4. Reverting to the original Eshelby 
problem, this amplitude can be interpreted as the product of a uniform plastic strain within a circular 
inclusion times its volume (or area in our two-dimensional framework). But the key point, is that the 
external stress field is independent of the specific shape of the inclusion and details of the re-arrangement 
- such dependencies only appear in higher multipolar orders. Nor does it make sense to distinguish 
separately volume and local plastic strain. In the sequel, conventionally the volume of the STZ can thus 
be chosen to be equal to the volume of the discretization scale, namely unity, and hence the amplitude r\ 
is a measure of the equivalent local plastic strain, denoted as "local slip" in the following. 

An elementary local transformation of amplitude r\ occurring at point Xq induces a local increase in 
the plastic strain field 

Ae p (x) = r)6(x - x ) (2) 

where 8 is the Dirac distribution. The residual stress field a e i is correspondingly modified by 

Aa e i(x) = r)G(x- x Q ) (3) 



This residual stress field itself encodes a memory of the set of transformations experienced by the medium, 
and long-range spatial correlation will naturally accumulate in the stress field. Thus at any stage of loading 
the local stress consists of a homogeneous macroscopic stress, S, added to the residual stress field, 

a(x) =E + a e i(x) (4) 

As local transformations can occur anywhere in space, each point x is attributed a threshold stress a 1 {x). 
Starting from rest, at each instant of time, we deduce the level of external loading needed to trigger a 
local transformation as 

E c = min[cr 7 (a;) - a e i(x)] (5) 

X 

The intrinsic local disorder of amorphous media implies that the slip amplitudes ry are randomly dis- 
tributed. Moreover, because our description is sought at a large scale no spatial nor temporal correlations 
will be considered. To account for the local reorganizations, it thus suffices to give the probability distri- 
bution function of their (scalar) amplitude. The only physical constraint is that they should be bounded. 
Similarly to the yield thresholds, the slip value r\ is drawn randomly from the uniform distribution, [Q,d] 
if not otherwise stated, d is a parameter of the model whose role will be discussed below. 

3.3. Extremal dynamics 

Quasi-static driving conditions are considered. Thus stress redistribution is considered to occur instan- 
taneously (no viscous effects are considered). However, one should finally specify the way the system 
is driven. A constant external stress is inappropriate as one expects the stress strain law to tend to a 
constant plateau stress. An imposed total strain is more relevant, but it has the drawback of possibly 
triggering multiple transformation zones. We have chosen the so-called "extremal dynamics", which con- 
sists of adjusting at each instant of time the external loading S so that one and only one transformation 
zone x*(t) is activated at each time t according to Eq. 5. This may involve a reduction in the applied 
macroscopic stress and strain. A unique temporal sequence of event is selected, from which it is easy 
to reconstruct any other driving mode assuming the same sequence of events. In this way, "avalanches" 
can be defined without ambiguity. Note however, that "time", t, is used here as a simple way of counting 
and ordering events. On average, time is simply proportional to the total plastic strain imposed on the 
system, (e p ) = tdL~ 2 /2. The plastic strain field is thus simply 

/ 
e p (x,t) = Y,r](t)6(x-x*(t)) (6) 

1 

This specific driving mode is the one that would result from an over-dampened viscous dynamics and 
infinitesimal strain rate driving. 

3.4. Stress redistribution function 

The Green function G can easily be computed to be G = cos(4#)/r 2 in polar coordinates for an infinite 
isotropic medium. The quadrupolar elastic kernel is obviously a key ingredient of the present model. As 
already discussed by Picard et al [52], the numerical implementation of this kernel is rather delicate. In 
order to deal with the bi-periodic boundary conditions we chose here to discretize first in the reciprocal 
space G(k,u>) = — cos(4u;) — 1, where k and w are the polar coordinates in the Fourier space, (Fig. la) 
then to Fourier transform to get the periodic kernel G(r,6) in the direct space (Fig. lb). In the present 
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Figure 1. Map of the quadrupolar elastic interaction used in the model. The discretization is performed in the Fourier space 
G(k,uj) = cos(4lj) — 1, where k and u> are the polar coordinates in the Fourier space (k x ,k y ) (a) and a subsequent Fourier 
transform gives in the direct space the Green function G(x, y) (b) satisfying the bi-periodic boundary conditions of the 
problem, (c) Spectrum of the eigenvalues of the elastic propagator. The horizontal tangent at the maximum eigenvalue zero 
can be associated to an abundance of soft modes (with null or close to zero eigenvalues). The associated eigen-modes are 
aligned along directions equal or close to ±45° . 

case, because of the bi-periodic boundary conditions, the elastic propagator is a circular matrix. The 
eigenvalues of this (convolution) operator are directly obtained by Fourier transforming the associated 
Green function. At the same time as the elastic propagator we thus obtain the spectrum of its eigenvalues, 
as shown in Fig. lc, without performing any further computation. 

The residual stress field resulting from an arbitrary cumulative plastic strain is easily written in Fourier 
space as 



a e i 



= G.e, 



(T) 



Hence integration over a long interval of time shows that a stationary state can be reached where a e i 
remains constant, while plastic strain accumulates over the modes associated with a zero eigenvalue for 
G. This shows trivially that uniform shear along lines at ±7r/4 orientation with respect to the x-axis, but 
with arbitrary values from line to line meets this stationary condition. However, as plasticity develops 
in our model by discrete events of random amplitude, it is not possible to fulfill such condition at each 
instant of time. The "granularity'' of the elementary events imposes a locality of plastic strain in real 
space, which in turn is non-local in Fourier space, and hence, stress fluctuations will occur for which it 
is essential to account faithfully for the vicinity of the spectrum of eigenvalues close to 0. This point is 
emphasized here because slight variations in the discretization of the Green function G which may appear 
as innocuous have a drastic influence in the long range scaling properties of the model. 

3.5. Slip amplitude distribution 



Let us come back to the meaning of the parameter d which is the only dimensionless parameter of the 
model. It quantifies the elementary plastic strain required to change the conformation of the local shear 
transformation zone, at the scale, a, of our discrete mesh in such a way that local yield stress correlation 
over time can be ignored. After one elementary local transformation event, the stress fluctuation for a 
site neighbor to the slip is of order d. This quantity is to be compared to the amplitude of the yield stress 
fluctuations, ?i. Hence the important dimensionless parameter is d/<,\. In the sequel, the yield stress 
amplitude ft is conventionally set to unity. Therefore, d remains the only free parameter of the model. 
For most of the simulations shown below, a value of d — 0.01 is chosen. 



In the following we will focus on the relation between the local plasticity events and the macroscopic 
plastic flow. We will show in particular that a scale invariant picture naturally emerges (and whose 
justification lies in the underlying criticality of the depinning transition[51]). It is therefore extremely 
important to identify those statistical features in order to be able to up-scale the description, and capture 
the transition toward a deterministic behavior at a large enough scale. 



4. Stress-strain relation 

As discussed above, the choice of an extremal dynamics corresponds to a quasistatic driving at a 
vanishing strain rate. While the strain rate is kept constant, the external stress is a highly fluctuating 
quantity. Far from being restrictive, this choice allows us to visit all successive configurations of the system 
and to reconstruct the response to any loading. The case of the response to a monotonically increasing 
stress is recovered by computing the maximum external stress over the past. However as this quantity is 
controlled by one single event in the past, it is a fragile quantity. In the following, the time evolution is 
cut into intervals corresponding to a number of event equal to the number of sites in the system (so that 
the mean plastic strain increment is d/2). The macroscopic stress is computed as the maximum external 
stress over each interval. Fig. 2 shows the evolution of this quantity, denoted here (£) as a function of 
the plastic strain for different values of the slip increment parameter d. 

It is to be noted that the stress-strain curve appears here to be a well-behaved quantity which is 
independent of the system size in the explored range (L = 16 to L — 256 for different values of the slip 
increment parameter d = 0.01 , d = 0.1 , d = 1 in Figure 2). Saturation occurs at a d-dependent value of 
the plateau stress. However, as can be seen in Figure 2, for higher values of d, a systematic size effect can 
be seen. It affects only the smallest sizes for d = 0.1, but appears to be very significant for sizes ranging 
from 16 to 256 when d = 1. 

A simple argument allowing us to rationalize the dependence of the plateau stress on the slip increment 
parameter d consists of regarding the elastic stress fluctuations as a mechanical noise comparable to 
a thermal one. Using this analogy early proposed by Sollich et al. [55] and more recently discussed by 
Lemaitre and Caroli[54], each local slip induces a stress fluctuation whose amplitude is proportional to d, 
which can be seen accordingly as an effective mechanical temperature. Following this interpretation, the 
higher d, the higher the mechanical noise and the lower the external stress needed to destabilize a shear 
zone. A clear limit of this interpretation is the high degree of spatial correlation of the elastic response to 
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Figure 2. Hardening curves: Averaged stress (£) vs. plastic strain e obtained for three values of the d parameter, (a! = 0.01 
(a), d = 0.1 (b) and d = 1 (c)). The curvature signs a transient hardening effect, i.e. the elastic limit (yield stress) increases 
with plastic strain. After the hardening stage, the stress saturates at a plateau value, the macroscopic yield stress Ey. 
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a local slip. Without pursuing further this interpretation, let us simply note that the systematic decrease 
of the plateau stress with the slip increment d observed in Fig. 2 follows qualitatively this picture. 

In Fig. 3a we show hardening curves obtained for a system of size L = 256 with values of the slip 
increment parameter ranging from d — 0.01 to d = 1. We then test in Fig. 3b a simple rescaling using 
the reduced variables (S)/Ey and (e p )/eo.5 where T,y is the plateau stress value and £0.5 is the plastic 
strain obtained at £ — 0.5£y. While not perfect the superimposition obtained with this procedure is 
reasonable and appears to gain in quality for small values of d. In Fig. 3c we show the scaling behavior 
of the plateau stress ay and the typical plastic strain £0.5- The latter shows a sublinear dependence to d, 
£0.5 oc d 085 while the former approaches unity, the upper limit of the threshold distribution with a close 
to square root dependence 1 — Sy oc <i ' 48 . 

A naive view of this behavior (leading to a square root scaling for both quantities) consists of noting that 
small values of d allows for a more precise exploration of the valleys of the disordered landscape. Multiple 
slips at the location are thus necessary before plastic activity jumps to another site. Since any time, a 
slip occurs, the internal stress is incremented by an elastic response of order d, the typical deformation 
can be seen as a first return of a biased random walk of slope d. 

The initial curvature of the stress/strain response is characteristic of a transient hardening behavior, 
also called micro-plasticity i.e. the progressive increase of the yield stress upon deformation. Such a phe- 
nomenon, traditionally attributed to dislocation entanglement or pinning by impurities in metal plasticity 
was first believed to be absent in amorphous plasticity [5]. Yet, the phenomenon of progressive densification 
of silicate glasses upon high pressure cycling can be interpreted as density hardening effect [64]. Moreover 
a clear effect of strain hardening upon shear cycling has been reported recently in a metallic glass [65]. 
Within the present model, it is possible to give a simple statistical interpretation of strain hardening 
phenomenon in absence of dislocations. As early discussed in Ref. [11] the increase of the yield stress can 
be related to the progressive exhaustion of the weakest spots of the materials. 

This phenomenon is clearly visible in Fig. 4 where we represented the evolution under shear of the 
distribution of local thresholds -P(c 7 ) and of the distribution of residual elastic stress P(a e i). While 
originally uniform in the interval [0, 1], the distribution of plastic thresholds is progressively shifted toward 
high values by exhaustion of the weakest sites before reaching a limit steady state distribution. In the 
transient regime, a self-organized-critical-like dynamics seems at work: weak sites are replaced by normal 
ones, thus inducing a systematic bias to the distribution. While less remarkable a similar trend may be 
extracted from the evolution of the stress distribution. A clear asymmetry is visible: large positive stress 
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Figure 3. (a) Averaged stress (E) vs. plastic strain e obtained for a system of size L = 256 with values of the slip increment 
parameter ranging from d = 0.01 to d = 1. The plateau stress value increases while the parameter d decreases, (b) same 
data after rescaling by the plateau stress value Ey and the plastic strain eo.5 obtained at (£} = 0.5Sy-. (c) Evolution of 
the plateau stress Eyand the typical plastic strain eo.5 with the slip increment d. 
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Figure 4. Evolution of (a) the distribution of microscopic plastic thresholds and (b) the distribution of internal stress under 
increasing plastic strain (from e p = 1/16 to e p = 8). Data are obtained with a system of size L = 64 with a slip increment 
parameter d = 0.1, averaged over 250 realizations. The horizontal black dot-dashed line gives the distribution of thresholds in 
the initial state (e p = 0). Due to the exhaustion of weak sites during deformation, the distribution of thresholds progressively 
shifts toward high values when e p inncreases and finally converges toward a limit distribution of typical width d. The internal 
stress distribution (here represented in semi-log coordinates) converges toward a Gaussian-like distribution with however a 
notable asymmetry: a fatter tail for the negative stress. 

values which favor local reorganization are less frequent than their negative counterparts. 



5. A non-persistent localizing behavior 

By the very definition of the model, the mean plastic strain increases in proportion to the number of 
events, or "time". As the mean plastic strain increases, the macroscopic stress reaches a plateau in the 
steady state, (dcr)/de = 0. Using classical results from continuum mechanics, for homogeneous media, the 
system is at the limit between a hardening behavior da/de > 0, expected to give rise to a homogeneous 
plastic strain, and a softening regime da/de < where strain should localize. A simple way to illustrate 
this property, is to note that a uniform (but of arbitrary magnitude) slip along a line oriented at ±7r/4 
with respect to the principal axes, will not induce any stress in the medium. Introducing n\ = l/y2(l, 1) 
and «2 = l/v2(l, — 1), and arbitrary scalar functions / and g, 



£ p (x) = /(as.ni) + g{x.n 2 ) 



(8) 



is a field which does not involve any elastic stress within the medium. In our model, the residual elas- 
tic stress is the only way to endow the system with memory, and hence, no limiting (nor amplifying) 
mechanism will act against (with) such degrees of freedom. 

5.1. Fluctuations of plastic strain and stress 



The presence of disorder gives rise to stress fluctuations around the average, which may affect this 
simple picture, and hence it is of interest to characterize the inhomogeneity of the plastic strain field. A 
first simple estimator is the variance of the entire plastic strain field, 



5e 2 p = [{el) - (e p f 



(9) 
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which is to be characterized as a function of the mean plastic strain (s p ). Indeed, on the one hand, for a 
plastic strain field which does not display any long range correlation, del, should saturate to a constant, 
independent of (e p ) and L. If, on the other hand, the plastic strain is localized on a simple slip system (a 
straight line spanning the entire system), then along this localization band, namely for L sites, the plastic 
strain amounts to E\ = L(e v ), whereas it is null over the rest of the system. In this case, the variance 
amounts to 5e p oc L(e p ) 2 . Therefore the scaling of 6e p with L and (e p ) is informative on the more or less 
uniform distribution of plastic strain. 

The evolution of the variance of the plastic strain as a function of the mean strain is shown in Fig. 5a in 
a log-log plot for systems of size L = 128 and various values of the slip increment parameter ranging from 
d = 0.01 to d = 1. A reasonable collapse onto a master curve can be obtained when using the reduced 
plastic strain e p /d}' 2 . As discussed above the quality of the rescaling procedure increases when the value 
of d decreases. 

A transient power- law like regime can be identified. An indicative straight line of slope 1.4 is shown in 
the transient regime. Then the strain fluctuation seems to slowly transit toward a diffusive like regime of 
slope close to unity. Again an indicative straight line of slope 0.95 is shown in this regime. Coming back 
to the above discussion we thus obtain a power law evolution of the variance of strain Se 2 oc (s p ) a with 
an exponent a m 1.4 in the transient regime and a ra 1 in a second diffusive-like regime. These results are 
thus intermediate between a homogeneous deformation regime (a = 0) and a full localization behavior 
(a = 2). 

In Fig. 5b we represent the evolution of the elastic stress variance. The contrast with plastic strain 
fluctuation is striking since we observe here a clear saturation after the transient regime. This means that 
plastic strain is dominated by stress free soft modes ±45°. These results help us to give a first picture of 
the behavior of the system. After spatial correlations have been established in a transient localized regime 
(high value of the exponent a), the system transits toward a diffusive regime where plastic deformation 
can be seen as a random succession of shear bands in the directions at ±45°. 
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Figure 5. (a) Log-log plot of the evolution of the rescaled variance of plastic strain 6Ep/d with rescaled plastic strain {e p )/d ' 5 
(L = 128, d = 0.01, 0.03, 0.1, 0.3, 1). A blue (or light-gray) straight line of slope 1.4 is shown as a reference in the transient 
regime, and a black straight line of slope 0.95 is shown in the "stationary" regime, (b) Evolution of the rescaled residual stress 
standard deviation 5cr e i/d -® 5 in the same coordinates. In contrast with strain, a real saturation regime can be identified. 
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5.2. Maps of plastic strain field 

Beyond this macroscopic behavior, a closer look at the local scale gives evidence for the development 
of strong spatial correlations. Fig. 6 presents a series of maps of plastic activity at different levels of mean 
plastic strain. Every map corresponds to the plastic strain accumulated during a finite strain window 
Ae p = 0.01. While the plastic activity is initially homogeneously distributed in space, a progressive 
development of correlations along the n\ and rt2 directions can be observed. Simultaneously, it appears 
that larger and larger regions remain still, i.e. without any plastic flow. As earlier noted, plastic strain 
appears as streaks aligned with either rt\ or ri2- Let us note that these patterns observed at the largest 
scale generate little stresses as they can grossly be described by the modes described in Eq. 8. 

Note that this localization appears not to be persistent. While the incremental plastic strain displays 
similar patterns at different times, shear bands are not superimposed, but they rather move through the 
system with statistically similar features past the initial transient at different deformation levels; they 
seem to diffuse in the system. The two left panels of Fig. 7 where we represented the plastic activity 
at e p = 0.95 and e p = 1.0 give an illustration of this competition between localization and diffusion. In 
the second map, we recover hints of the first one but most of the plastic activity has moved elsewhere. 
The last panel of Fig. 7 proposes a striking comparison with recent numerical results by Maloney and 
Robbins[23,26]. These authors study the behavior of a two-dimensional Lennard- Jones glass under uni- 
axial compression, and periodic boundary condition for the stress and strain field, a situation which is 
quite comparable to our own boundary conditions. They show the vorticity of the displacement field, in 
order to highlight the zones of concentrated plastic strain. With this representation, a slip line along the 
ri\ (resp. n^) direction appears as a line of concentrated negative (resp. positive) vorticity, whereas the 
local density of plastic strain does not distinguish between both directions. 

The growth of the correlation length takes place in the initial transient up to the stage where it reaches 
the system size, and steady state is established. Up to now, we essentially focused on this initial stage. 
However, since we observed the formation of large scale structures, it is of interest to characterize the 
signature of the largest structures in the steady state. In the following sections space correlations in the 
stress and plastic strain fields are discussed. 



6. Anisotropic plastic strain correlations 

We now give a quantitative analysis of the anisotropic plastic strain correlations. A first method consists 
in considering the spatial distribution of successive plastic events P(Ax, Ay). A map of this distribution 
is presented in Fig. 8(a) for the upper right quadrant. Without surprise, we obtain a panache that obeys 
the quadrupolar symmetry of the elastic kernel: a plastic event is far more probable in a direction at ±7r/4 
degree than along the principal axis. The width of this panache can be quantified. Fig. 8(b) shows the 
standard deviation sd& y of the distance Ay for a given value of the distance Ax between two successive 
plastic events. We obtain a nice scaling relation sd^y <x Ax 1 * where ( ps 0.65. Note that a similar value 
was obtained in the anti-plane shear case discussed in Ref. [11]. Such a value corresponds also to the 
roughness exponent of minimal path in a random potential [66] which can easily be proven to be relevant 
for anti-plane plasticity with quenched random yield stress. The fact that the value of this exponent is 
less than unity indicates an important property of the plastic strain patterns: the aspect ratio decreases 
as Ax^ _1 so that asymptotically, shear is concentrated on straight bands of vanishing relative width. 

In order to characterize the correlations of the plastic strain field, the most natural analysis is to 
compute the power spectrum of the plastic strain field. In Fig. 9a, an ensemble average of 2D power 
spectra of the strain field is shown. This map, here encoded in logarithmic scale, exhibits a marked 
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Figure 6. Maps of cumulated plastic activity during a finite strain window Ae p = 0.01 corresponding to increasing levels of 
deformation: e p = 0.01,0.10,0.20,0.30,0.50,1.00. The system size is L = 256. While random at the beginning, the plastic 
activity progressively localizes along n\ and ri2 directions. 
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Figure 7. (a) Map of cumulated plastic activity during a finite strain window Ae p = 0.01 taken at e = 0.95 with L = 256, 
d = 0.01. A clear localization of the plastic deformation is observed, (b) same as (a) but taken "later" at e p = 1.0. Again 
plasticity localizes along elongated structures at ±45° but the patterns are markedly different: localization is not persistent, 
(c) for comparison, reproduction of a strikingly similar map of plastic activity (vorticity of the displacement field) recently 
obtained by Maloney and Robbins[26]) on a 2D Lennard- Jones glass under compression. 
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Figure 8. (a) Map of P(Ax, Ay), the probability that two successive plastic events are separated by distances Ax and Ay 
along x— and y-axis for a system of size L = 256. Preferential angles at ±7r/4 are clearly visible, (b) Standard deviation 
sd&y of the separation of two consecutive plastic events along the y— axis for a given value of the separation distance Ax. 
A power law behavior sd^ y <x Ax<> is obtained with f Ri 0.65. The fact that f < 1 shows that the relative width of those 
shear bands tends to for large system sizes. 

fourfold symmetry, with obvious preferential directions along n\ and rii- To be more quantitative, we 
now consider the scaling of power spectra cuts along the preferred directions and along the principal axes, 
as shown in Fig. 9b. Both of these spectra can be characterized by a power-law behavior with exponents 
a T /4 ss 1.7, ao ~ 0.3 where the indices denote the orientation of the wave vector with respect to the x- 
axis. Such behaviors can equivalently be characterized as self-affine, with respective roughness exponents 
Qk/a ~ 0.35 and Co ~ —0.35. It is however to be emphasized that the most important feature of power 
spectra along directions 8 ^ ±7r/4, is that their magnitude is orders of magnitude lower than for 8 ± n/4. 

In the case of anti-plane shear, a previous study [11] also revealed an anisotropic plastic strain field in 
the steady state. In that case, the power spectra of e p showed a different scaling with the wave-number 
parallel or perpendicular to the orientation of the shear band. 

In their recent study of a Lennard- Jones glass under compression [23,26], Maloney and Robbins dis- 
cussed the anisotropic scaling of plastic strain. As shown in Fig. 7, when representing the vorticity of 
the displacement field, they obtained maps of plastic activity strikingly similar to ours. Looking at the 
correlation of the vorticity they obtained a power spectrum with a four-fold symmetry: 



S(q, 8) ex a(8)q~ a(e) , a{8) = 0.68 - 0.5 cos(40) 



(10) 



and thus in particular q t /4 = 1.18 and ao = 0.18, to compare to the above reported values. In the latter 
expression, a(8) is unspecified. 

Considering the power spectra of the cumulated plastic strain, we expect to observe a directly com- 
parable quantity. Indeed, we do observe a direction dependent scaling exponent. In Fig. 9 we reported 
the scaling behavior obtained along cuts in different directions while the angular dependence of S(q, 8) is 
shown in Fig. 10. Again most of the plastic activity appears to be concentrated along the directions at 
±7r/4. Except in the close vicinity of that preferred direction, the value of the scaling exponent remains 
close to ao, and only shows a slow increase when the direction of the cut approaches the diagonal. 

Representing part of the same results as a function of the polar angle, we observe that an equation 
comparable to the previous, gives a good account of our data 
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Figure 9. (a) Power spectrum of the plastic strain using a decimal log scale. A clear quadrupolar symmetry is obtained, (b) 
Cuts along the axis at 8 = and 8 = 7r/4. Different power law behaviors S(q) oc q~ a are obtained between the extreme 
cases ay ~ 0.33 and a^u £S 1.67. L = 256 and d = 0.01. 



S(q, 9)= A (q/q r a{e) , a(9) = a n/4 - K /4 - «o)| cos(20)| 



0.4 



(11) 



where we considered the case d = 0.01. The angular dependency of the power spectrum is written in 
a slightly different manner, but both share the same fourfold symmetry. At a fixed wavevector modu- 
lus, Eq. 11 proposes an amplitude of Aq$ whereas Maloney and Robbins propose A(9). However the 
difference is marginal, as the lower exponents correspond to much smaller amplitudes. 
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Figure 10. Angular dependence of the power spectrum obtained for two values of the wavenumber: q = 16 (a) and q = 64 
(b) with the parameters L = 256, d = 0.01. The symbols represent the numerical data; The dashed lines corresponds to the 
best fit obtained with a scaling form s(q,8) ex A(q/qo)~ "' ' with an exponent ce(8) = 1.67— 1.34 cos(20)| 0,4 exhibiting a 
quadrupolar symmetry. 
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7. Conclusions 

The plasticity of amorphous media was addressed through a meso-scale modeling in order to highlight 
the large scale mechanical behavior. This modeling was constructed in order to retain the most impor- 
tant features of a macroscopic description and compared favorably with much richer descriptions. The 
macroscopic behavior tends to approach an elastic/perfectly plastic law. The progressive approach to a 
constant yield stress is interpreted as a statistical selection of high threshold local configuration. In the 
perfect plasticity regime, strain tends to localize along the direction of maximum shear as could have been 
anticipated. However, the fact that elementary plastic events are localized in space give rise to specific 
features in the spatial correlation of plastic strains (with self-similar power spectra) , and fluctuations of 
macroscopic yield stress occurring at large scales. Let us finally note that this model has recently been 
shown to exhibit a non-trivial avalanche behavior [67]; moreover the influence of the initial state of the 
medium has been shown to favor an early localization of the strain onto a unique and persistent shear 
band[68]. 

An extension accounting for the coupling between deviatoric and volumetric strain remains to be 
investigated. It would also allow to study the hardening behavior while paying attention to the potential 
apparition of an anisotropic texture as recently discussed in [39] . 

We are thankful to C.E. Maloney and M.O. Robbins for the authorization of replicating their figure in 
Fig. 7 from Ref [26]. 
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